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1  Introduction 


In  tins  work  wc  investigated  several  geometrical  issues  arising  in  high-frequency  electro¬ 
magnetic  (EM)  calculations  on  aircraft  modeled  with  nonuniform  rational  bi-cubic  splines 
(NURBS).  An  aircraft  model  represented  by  NURBS  consists  of  a  collection  of  surface  patches 
each  of  which  is  described  by  rational  bi-cubic  splines.  Such  representation  of  surfaces  is 
widely  accepted  in  design  and  manufacturing  industries  and  utilization  of  the  already  exist¬ 
ing  models  foi  electromagnetic  analysis  is  a  highly  desirable  alternative  to  the  very  expensive 
a.nd  time  consuming  task  of  building  anew  a  suitable  computational  mesh.  However,  one 
of  the  main  difficulties  in  modeling  and  computing  wave  propagation  over  aircraft  repre¬ 
sented  by  NURBS  is  connected  with  the  fact  that  while  each  surface  patch  has  a  smooth 
parametmation,  the  overall  smoothness  may  not  be  available  because  on  typical  models  the 
patches  do  not  match  along  their  boundaries  within  the  accuracy  limits  required  for  the  EM 
analysis.  In  other  words,  such  models  have  artificial  discontinuities.  Obviously,  if  »  model 
is  to  be  used  for  the  EM  analysis  these  discontinuities  should  he  recognized  and  accounted 
for.  Analytically,  these  discontinuities  result  in  discontinuities  in  coefficients  of  differential 
equations  describing  the  wave  propagation. 

A  possible  way  to  deal  with  this  problem  is  to  perform  a  pre-processing  step  on 
which  a  given  model  is  modified  so  that  the  rather  stringent  geometric  requirements  arc* 
met.  Unfortunately,  very  few  techniques  for  performing  such  a  step  are  currently  available 
and  in  many  cases  an  extensive  manual  and  time-consuming  effort,  for  '‘healing’'  a  model  is 
required. 

In  this  work  we  investigated  an  alternative  strategy  for  dealing  with  this  problem. 
In  our  approach,  rather  than  performing  an  expensive  general  pre-processing  step,  wc  are 
recognizing  such  discontinuities  within  the  LM  analysis  process  and  provide  means  for  dealing 
with  them.  Specifically,  for  the  particular  problem  of  surface  ray-tracing  required  for  surface 
diffraction  analysis  we  developed  several  techniques  for  continuing  wave  propagation  paths 
(geodesics  in  the  surface  metrics)  across  the  boundaries  of  patches.  These  techniques  are 
based  on  application  and  extension  of  geometric  results  on  manifolds  with  singular  metrics 
[2].  [3],  [4].  An  algorithm  realizing  this  approach  to  rav-tracing  lias  been  also  developed  and 
implemented  in  a  system  of  computer  codes.  The  codes  have  been  tested  on  some  simple 
surfaces  and  on  the  fighter  VFY  21$  model  in  the  ICES  12$  format  (sec  section  2  for  details 
on  this  format).  The  results  of  the  tests  confirmed  the  validity  of  our  approach. 


Future  work  may  be  an  expansion  of  the  developed  codes  to  include  more  complex 
propagation  paths  containing  surface  and  spatial  segments  as  well  as  edge  and  corner  diffrac¬ 
tions.  reflections,  etc.  These  developments,  in  conjunction  with  an  implementation  of  the 
Geometric  and  Uniform  Theories  of  Diffraction  (GTD  and  UTD,  respectively)  based  on  the 
developed  geometric  computational  techniques,  should  lead  to  a  new  set,  of  computational 
codes  providing  better  accuracy  of  high-frequency  KM  analysis  titan  that  available  wit  h  the 
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currently  existing  codes. 

The  theoretical  results  obtained  under  this  effort  are  currently  being  prepared  for 
publication. 


2  IGES  Format  and  Surfaces  Represented  by  NURBS 

A  representation  of  surface  data  with  NURBS  can  be  organized  in  many  different  ways. 
In  this  work  we  utilized  the  representation  provided  by  the  so-called  Initial  Graphics  Ex¬ 
change  Specification  (IGES)  format,  version  5.3  [!].  The  IGES  format  is  a.  large  collection 
of  information  structures  for  digital  representation  and  exchange  of  product  definition  data, 
including  curves  and  surfaces.  This  format  evolved  over  a  period  of  almost, 20  years  and 
it  is  accredited  as  an  American  National  Standards  Institute  (ANSI)  Standard.  Because 
of  the  strong  interest  in  utilizing  aircraft  models  in  IGES  formal  by  the  Computational 
Electromagnetics  community,  it  is  important  to  investigate  the  issues  of  implementing  wave 
propagation  techniques  on  such  models. 


in  the  IGES  format  different  data  structures  for  describing  the  geometry  of  objects 
are  distinguished  by  the  “Entity  Type  Number”.  The  specific  IGES  entity  used  in  this  work 
for  aircraft  model  representation  is  Die  IGES  128.  A  complete  description  of  all  entities  in 
the  IGES  can  be  found  in  [l] . 

In  the  IGES  128  format  a  surface  is  represented  as  an  indexed  collection  of  surface 
pieces  each  of  which  is  a  grid  of  surface  patches  constructed  from  rational  functions.  The 
locations  of  the  nodes  of  the  gricl(s),  the  degrees  of  the  polynomials  defining  the  surface 
patches,  and  their  corresponding  coefficients  arc  stored  in  the  appropriate  data  structure. 
The  format  of  the  entity  follows  specific  rules  that  must  be  used  by  an  application  in  order 
to  retrieve  from  a  given  data  structure  the  required  geometric  and  analytic  information. 

A  complete  description  of  a  surface  is  provided  by  the  parametric  equations  of  a 
surface,  that,  is,  the  coordinate  functions  describing,  the  surface  a.s  a  submanifold  in  If.  In 
our  case,  the  parametric  equations  of  each  surface  patch  are  recovered  from  the  IGES  data 
structure  describing  the  surface. 
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3  Geodesics  on  Surfaces  Represented  by  NURBS 


High  frequency  approximations'  require  often  the  knowledge  of  surface  geodesics  along  which 
surface  waves  propagate;  see  Kouyoumjiau-Pathak-Burnside  [5],  Pat hak- Wang  [G]  [7],  [$} 
and  other  references  there.  For  surfaces  represented  by  NURBS,  before  one  can  proceed 
with  the  construction  of  surface  geodesics,  it  is  necessary  to  determine  the  coefficients  of 
the  differential  equations  defining  the  geodesics.  These  are  constructed  from  the  parametric 
equations  of  the  surface  as  follows. 

Let  5  be  a  surface  in  R?  given  as  a  collection  of  patches.  Each  patch  is  parametrized 
by  coordinates  (u,v)  which  arc  the  Cartesian  coordinates  in  a  domain  D  on  a  it, v  -  plane. 
On  each  patch  the  surface  is  represented  bv 

r(v,v)  =  (*(v,v),y(«,v),x(tt,t>))t  (u,v)  €  D,  (ij 

where  the  functions  x(u.  v)..y(»,v)  and  s(ti,v)  are  rational  functions  determined  by  the  ICF.S 
data  structure.  The  first  fundamental  form  (that  is,  the  metric)  of  S  is  a  quadratic  form 

ds*  —  E(n,v)du2  +  2F(u,  v)dudv  4-  C(u,v)dv 2, 


where 

E  —  QL  —  p  —  —  C  — 

du  du  ()u  c)v '  <)v  c)v 

Using  the  coefficients  E,  F  and  C  the  Christoffe]  symbols  of  the  second  kind  are  constructed 
according  to  the  formulas; 


and 


1  u  — 


r1  — 
1 22  — 


GEU  -  2  FFV  +  FE, 
2(EG-  F2) 

pi  _  GEv  —  FGU 

12  _  2 {eg -ny 

2GFV  -  GGu  -  FGV 
2 (EG  -  E2) 


r?i« 

2  EFU  -  EEV  -  FEU 

2  [EG-E2] 

(2) 

r?a  = 

GEU  -  FEV 
'  2 (EG-F2)’ 

(3) 

__  EGV  -  2  FFV  +  EG,, 

(4) 

^22 

2  (EG-F*) 

r2 

1  21 

-  r2 

-  1 

where  Ev,...  denote  the  corresponding  partial  derivatives. 


The  system  of  ordinary  differential  equations  defining  a  geodesic  is  given  by 


v>  =  />, 

p>  =  -V[y  -  2i"uPq  - 


vt  =  <7. 

¥  =  “I  -  2VUpq  -  l'2, ci2. 


(5) 

(0) 


Observe,  that  the  coefficients  of  the  system  (5),  (C>)  are  computed* only  for  a  specific 
patch  and  will  change  when  we  move  to  a  different  patch.  Obviously,  in  order  for  the 
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Ghristoflel  symbols  to  remain  continuous  across  the  boundaries  of  adjacent  patches,  these 
patches  must  match  up  to  the  second  order  of  smoothness.  As  it.  has  been  explained  in  the 
introduction,  quite  often,  this  condition  is  not  satisfied  and,  consequently,  the  right  hand 
side  of  (5)  and  the  coefficients  in  (6)  are  discontinuous.  We  will  return  to  this  issue  in  the 
next  section,  after  we  complete  the  discussion  of  the  initial  conditions  for  the  svstern  (5) 
(6).  - 

Different  types  of  boundary  value  problems  can  be  considered  for  the  system  (5).  (0). 
For  a  point  source  of  radiation  on  a  surface  it  is  natural  to  consider  the  initial  boundary- 
value  problem  (IVP)  in  which  the  location  of  the  source  and  the  direction  of  propagation  are 
specified.  Thus,  the  user  must  specify  the  initial  point  P  on  the  surface  (the  wave  launching 
point  where  the  geodesic  should  start)  and  the  initial  direction  T  (of  unit  length)  in  which 
the  geodesic  should  go. 


The  initial  conditions  require  some  explanations.  The  first  difficulty  is  that,  typically, 
the  initial  point  is  given  as  a  point  in  723  (even  if  it  is  on  the  surface),  but  tin;  system  (5),  (fi) 
is  in  terms  of  the  local  coordinates  u,  v.  Therefore,  for  a  given  initial  point,  say,  A'o,  Vy./?0, 
one  needs  to  determine 

(a)  the  surface  patch  it  belongs  to  and 

(b)  the  local  coordinates  «0,t'o  such  that  -Y{u0,w0)  =  A'0,  F(itu,  v0)  =  Y0.Z{vv.v0)  =  ZQ. 

The  step  (a)  is  accomplished  by  a  search  algorithm  and  the  step  (b)  by  a.  gradient 
method. 


Finally,  we  need  to  set  up  the  initial  direction  in  terms  of  local  coordinates.  Suppose 
the  initial  direction  is  given  as  a  vector  T  €  If?.  Once  the  point  (u0.v0)  is  known,  the 
derivatives  r„(«o,vo)  and  r„(f/o,vo)  can  be  computed.  Then  the  vector  T  is  decomposed  in 
these  two  vectors  and  we  obtain  its  corresponding  component  in  required  form.  (To  make 
sure  that  T  is  tangent  to  the  surface,  we  take  as  T  only  the  projection  of  the  initially  given 
vector  onto  the  plane  spanned  by  r„(wo,vo)  and  r„(uo,vo).)  Further,  T  should  be  normalized 
to  have  unit  length.  The  corresponding  two  numbers  (that  is,  the  coefficients  po  and  </0  of 
the  expansion  of  T  in  ru(u0.  t'o)  and  r,,(«0, n0)  provide  the  remaining  initial  conditions  for 
the  system  (5),  (6),  that  is. 

p(wo.Vu)  =  py,  (7) 

*  '/(«IJn  '’o)  =  (]U.  (N) 


G 
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4  Construction  of  Numerical  Solutions  to  (5)  -  (8) 

1.  Let  P  be  the  patch  containing  the  initial  point  u0,z>0.  On  P  the  system  (5)  -  (8)  is  solved 
by  a  high  order  ordinary  differential  equations  solver.  The  solution  (that  is,  the  geodesic) 
is  continued  until  it  reaches  the  boundary  of  P.  Let  Q  be  the  last  point,  on  the  geodesic 
in  patch  P.  In  order  to  continue  the  geodesic  past  the  boundary  of  that  patch,  we  first, 
identify  the  patch  of  the.  surface  that  also  contains  Q  or,  if  such  a  patch  can  not  be  found, 
we  identify  the  patch  for  which  the  distance  to  Q  is  minimal  among  all  patches  (except  P). 
If  this  minimum  is  within  the  user  specified  error  tolerance  and  if. such  patch  is  unique  then 
we  have  a  patch  into  which  the  geodesic  is  to  be. continued.  If  there  are  several  patches 
on  which  the  same  minimal  distance  is  realized,  additional  steps  are  required  to  make  the 
proper  selection.  We  will  describe  these  steps  below  in  subsection  2.  If  no  suitable  patch  is 
identified  the  algorithm  stops. 

Suppose  that  the  algorithm  found  a  unique  patch  into  which  the  geodesic  should  he 
continued.  Denote  this  patch  by  P'  and  let,  Q'  be  the  point,  on  P'  closest  to  Q  (including 
the  case  when  Q  =  Q',)  The  semi-tangent  to  the  geodesic  at  the  point  Q  (on  the  side  of  P) 
defines  a  direction  T'  in  which  the  geodesic  should  be  continued.  If  V  is  tangent  to  P'  at.  Q' 
then  Q'  and  T'  are  taken  as  the  new  initial  conditions  in  patch  P' .  The  coefficients  of  the 
system  (5)  -  (6)  are  redefined  in  this  patch  and  a  new’  piece  of  the  geodesic  is  constructed. 


If,  however,  the  vector  T'  is  not  tangent,  to  P\  we  perform  the  following  operations, 
hirst,  we  check  that  the  tangents  to  the  boundary  curves  of  P  and  P'  passing  through  Q 
and  Q\  respectively,  are  parallel.  Tf  this  is  the  case,  then  the  vector  T’  with  the  base  point 
Q'  is  rotated  ill  the  plane  perpendicular  to  that  tangent  until  it  becomes  tangent  to  the 
patch  P' .  If  this  is  possible,  then  the  corresponding  new  direction  is  the  direction  in  which 
the  geodesic  is  continued.  This  procedure  implements  the  so-called  “equal  angle”  condition 
based  on  geometric  properties  of  geodesics  on  manifolds  with. singular  metrics  [2],  [9],  It  can 
be  shown  that  this  is  the  correct  continuation  of  a  geodesic  in  the  presence  of  discontinuity. 
This  procedure  allow’s  to  handle  also  the  case  when  the  discontinuity  is  due  to  the  presence 
of  an  actual  edge  (say,  an  edge  of  a  wing  on  an  aircraft). 

If  the  tangents  to  the  boundary  curves  of  P  and  P'  passing  through  Q  and  Q'. 
respectively,  are  not  parallel  (within  specified  error  tolerances)  the  algorithm  stops. 

2.  Let  us  now  return  lo  the  case  when  there  is  more  than  one  patch  at  the  same 
distance  from  Q.  In  this  case,  the  surface  has  either  a  self-intersection,  or  it.  touches  itself 
at  Q.  or  the  point.  Q  is  a  corner  point  for  several  patches.  The  eases  of  a.  self-intersection 
and  self-touching  are  not  typical  for  aircraft  surfaces  and,  for  that  reason,  not  considered.  In 
the  case  where  several  patches  have  Q  as  their  corner  point  the  following  procedure  is  used 
to  identify  the  proper  cont  inuation  of  the  geodesic.  First,  for  each  patch  with  the  vertex  Q 
we  compute  the  surface  angle  between  semi-tangents  emanating  from  Q  ami  tangent  to  the 
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boundary  curves  emanating  from  Q.  Let  K(Q)  be  the  sum  of  these  angles  taken  over  all 
patches  with  the  vertex  Q.  Let  l  be  the  vector  tangent  to  the  already  constructed  geodesic  in 
patch  P  and  directed  “into”  the  patch  P.  Define  a  vector  V  tangent  to  some  patch  adjacent 
to  Q  and  such  that  the  angle  between  t  and  T  measured  on  the  surface  consisting  of  patches 
adjacent  to  Q  is  equal  A'(Q)/2.  Such  a  vector  is  uniquely  defined  and  its  direction  is  taken 
as  the  direction  for  continuing  the  geodesic.  This  procedure  is  also  based  on  properties  of 
geodesics  on  manifolds  with  singular-  metrics  (2],  [3];  see  also  . 

In  the  beginning  of  the  algorithm  the  user  specifies  the  length  of  the  required  geodesic 
and  the  above  steps  are  repeated  until  either  the  required  length  is  attained  or  the  algo¬ 
rithm  stops  because  the  surface  discontinuities  do  not  allow  a.  reasonable  continuation  of  the 
geodesic. 


Figure  1:  A  propagation  path  on  the  VFY213  fighter  modeled  with  NURBS 
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